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Abstract: We study the dynamics of active-sterile neutrino oscillations in the early 
universe using full momentum-dependent quantum-kinetic equations. These equa¬ 
tions are too complicated to allow for an analytical treatment, and numerical so¬ 
lution is greatly complicated due to very pronounced and narrow structures in the 
momentum variable introduced by resonances. Here we introduce a novel dynamical 
discretization of the momentum variable which overcomes this problem. As a re¬ 
sult we can follow the evolution of neutrino ensemble accurately well into the stable 
growing phase. Our results confirm the existence of a “chaotic region” of mixing pa¬ 
rameters, for which the final sign of the asymmetry, and hence the SBBN prediction 
of '^He-abundance cannot be accurately determined. 
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1. Introduction 


Recent observations from atmospheric and solar neutrinos from Super Kamiokande 
(SK) [1] and Sudbury Neutrino Observatory (SNO) [2] prefer conventional, mainly 
active-active mixings between ordinary electron, muon and tan neutrinos as the so¬ 
lutions to the observed flux dehcits. However, while pure active-sterile solutions 
are disfavoured, sizable admixtures of sterile states in the observed fluxes are still al¬ 
lowed [3]. New sterile states are also needed in order to explain the LSND anomaly [4] 
by neutrino physics. Moreover, in a recent analysis of the so called “2-1-2”-models^ 
for four-neutrino mixing, the solutions with nonzero sterile neutrino components 
were found to provide best global hts to the solar, atmospheric and reactor data [6]. 
Corresponding limits on active-sterile mixings are quite generous, for example 

sin^ < 0.48 (Atmospheric) 

sin^^e. <0.72 (Solar LMA), (1.1) 

where LMA refers to the Large Mixing Angle solution for the solar neutrino dehcit. 

Active-sterile mixing, if realized, would have several interesting effects in astro- 
physical settings [7] and in particular for the evolution of the early universe [8-20]. 
In particular, sterile states could be brought into thermal equilibrium by mixing be¬ 
fore nucleosynthesis, so that the resulting anomalous increase in the expansion rate 
of the universe would lead to overproduction of helium in disagreement with the 
observations [8-13]. Excluding the parameter sets leading to equilibration provides 
useful bounds on active-sterile mixing. Indeed, from results of ref. [11] one can infer 
that the sterile components in mass eigenstates responsible for atmospheric anomaly 
and LMA are constrained by 

sin^ < 0.013 (Atmospheric) 

sin^ 9es < 0.026 (Solar LMA) (1.2) 

where we used a SBBN-limit of < 3.4 for the number of effective neutrino degrees 
of freedom [11,21]. These numbers correspond to the “light” case where the sterile 
component is the heavier of the mixing states; in the opposite, “dark” case, the 
corresponding limits are even stronger by one to two orders of magnitude. In either 
case the cosmological constraints are more stringent than those obtained in terrestrial 
laboratories. 

While quite generic, the cosmological bounds (1.2) do depend on some simple 
prior assumptions, the most important of which is that the primordial lepton asym¬ 
metry is not anomalously large [14]. It may be possible to circumvent them in more 
complicated mixing schemes involving more sterile states (for a recent discussion, 
see [22]). In a particular attempt it was shown in ref. [15] that for a large negative 

^Alternative “3-1-1” models fit less well with short baseline reactor experiments [5]. 
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and a small enough mixing angle (so that the equilibration bounds of [11] can 
be avoided), say on the Vr — z^s'-sector, resonant oscillations may trigger a rapid 
growth of the tan lepton asymmetry. This asymmetry could then become very large, 
Lj. ~ 0{1) and change the natural prior condition of no signihcant initial lepton 
asymmetry for a mixing in any other sectors. Indeed, if created early enough, such 
an asymmetry could suppress — z/^-oscillations [14,23], obviating the bounds (1.2) 
for z/^ — z/g-mixing. However, as this scenario requires that the sterile state is the 
lighter one (resonance condition), and that the mass splitting is large (to cre¬ 

ate L-r early enough), it is essentially excluded by the recent solar and atmospheric 
neutrinos combined with the direct constraint on the electron neutrino mass [24]. 

Even with the mechanism of ref. [15] by and large excluded, there can be impor¬ 
tant effects due to L-growth with smaller For example, a large homogeneous 

electron neutrino asymmetry would directly influence the nucleosynthesis pre¬ 
diction for the hehum-4 abundance through the reactions 

n -f e"*" ^ p + Ue 

p + e~ ^ n + Ue , (1.3) 

which keep the neutron-to-proton ratio in equilibrium at early times. Changing the 
electron neutrino abundance could tilt the balance of these reactions with the pos¬ 
sibility of either increasing Yue (negative and hence strengthening the bounds 
(1.2), or decreasing Yne (positive L^J, leading to weakening of the bounds (1.2). 

The physics involved with the resonant growth of the asymmetry is very compli¬ 
cated because of several vastly different physical scales both in the temporal direction 
and in momentum variable, which effect the evolution of the ensemble in an essential 
way. The relevant QKE’s are nonlinear and strongly coupled internally via the asym¬ 
metry term, and therefore no useful analytical approximation exists for the problem^. 
Numerical solution is also greatly complicated by the various different physical scales, 
and while the phenomenon of asymmetry growth is fairly well established by now, 
the determinacy of the hnal sign of the asymmetry is still not well understood. The 
ambiguity was hrst observed in [16], and in [17,19] it was shown that this “chaoticity” 
occurs in certain well dehned region of mixing parameters. References [16,17,19] em¬ 
ployed a numerical solution for the momentum averaged approximation of the QKE’s 
however, and their results were challenged by ref. [25], who claimed that the hnal 
asymmetry is always fully determined by the initial conditions. While the analysis 
of ref. [25] is now discredited [26], it still remains true that no numerically reliable 
momentum dependent analysis has so far shown the existence of the chaotic region.^ 

^Some confusion in the field was caused by a recent analytical treatment [25], whose results 
contradicted earlier numerical results [16-18]. Disagreement was eventually clarified in favour of 
numerical work in [26]. 

^Although ref. [27] found support for chaoticity, their conclusion is not firm because of the 
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It is important to settle the issue because, as mentioned above, large could sig¬ 
nificantly alter the helium-4 abundance. Indeed, if the sign of were found to be 
chaotic, then Yhc could not be reliably determined from BBN [17]. 

In this paper we study the dynamics of active-sterile neutrino oscillations in 
the early universe using full momentum-dependent quantum-kinetic equations (in 
homogeneous space). Our main technical improvement over the previous works is the 
introduction of a novel dynamical discretization method for the momentum variable 
which enables us to model accurately the density matrices over the entire momentum 
range, including the extremely pronounced and narrow structures close to resonant 
momenta. On physical context, our results confirm the existence of a ’’chaotic” region 
of mixing parameters, for which the final sign of the asymmetry is not deterministic. 
It also confirms the expectation [17,27], that the size of the chaotic region is somewhat 
smaller than indicated by the momentum averaged code [17]. These results can be 
seen as giving more strength on the cosmological constraints on neutrino mixing 
parameters. 

The paper is organized as follows. In section 2 we set up the quantum kinetic 
equations for the problem. In section 3 we introduce the novel dynamically adjusted 
discretization of the momentum variable and in section 4 we set up our kinetic 
equations with this parametrization. The method we develop allows us to have 
enough resolution to solve the problem with relatively small number of momentum 
bins (we use at most 400 bins). In section 5 we display our numerical results for 
the key variables driving the oscillation and discuss the numerical stability of our 
solutions. Finally, in section 6 we present our results on the asymmetry growth and 
oscillations as well as the interpretation of the physical consequences, and the section 
7 contains a summary and outlook. 


2. Quantum Kinetic Equations 

In the early Universe the oscillating neutrinos experience frequent scatterings which 
interrupt the coherent evolution of the state and introduce collisional mixing between 
different momentum states. The mathematical formalism which can incorporate all 
these features has been developed elsewhere [10,11,28-30], and it takes the form of 
the quantum kinetic equations for the reduced density matrices for the neutrino and 
antineutrino ensembles. We parametrize the density matrices by the Bloch vector 
presentation [11]: 


Pi/ = -/o(-Po + P ■ <^) ; pp = 2/o(-Po + P ■ <^); (2.1) 

reported loss of the numerical accuracy of their methods when approaching the potentially chaotic 
region. 
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where a are the Pauli spin matrices and /o = (1 + exp(p/T))“^ is the Fermi-Dirac 
distribution function without chemical potential. Then to the the lowest order in the 
interaction [10,28-30], the QKE’s for the density matrix take the form 

r 

dfPo [feq Paa\ 

Jo 

dtPx = -VzPy - DPx 
dtPy = V^Px - V^Pz - DPy 

dtPz = VxPy + -T [feq ~ Paa] : (2-2) 

Jo 

where dt = dt — Hpdp and H is the Hubble expansion factor. The rotation vector 
V has the following components {Vy = 0): in x-direction one has only the vacuum 
contribution 

dm? 

14 = ^—sin 26^0 , (2.3) 

2p 

whereas in the ;2-direction also the matter induced effective potential contributes: 

14 = Ho + Hi + Vl, (2.4) 


where 


Ho = 
H4 = 

Hi^’" = 

Hl = 


5m^ 

2p 


cos 26^0, 


-14y2^^iV..^p T{j + j cos^ Pw K, + n^J) 
7v^C(4) Gf 


2 C(3) M| 

V2GfN^lG\ 


N^p T 


(2.5) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 


Here is the photon equilibrium number density, is the normalized to unity 

(in equilibrium) actual neutrino (antineutrino) number density, and the effective 
neutrino asymmetries L*-") given by 

~ (7 4” ^ + (f ~ 2 sin^ Lp — + 2L^^ + (2.9) 

lG) = + (2.10) 

Lir) = lG) -L,-L,^ + L,^. (2.11) 


where Lf = {rif — nj)Nf/N^. 

The repopulation term — paa) in the diagonal part of (2.2) is written in the 
relaxation time approximation [31]. The distribution /g^ in the repopulation term is 
the usual equilibrium Fermi-Dirac distribution 


feq = 


1 -f CT 


J± 

T 


( 2 , 12 ) 
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and the reaction rates are 

r = C^GlpT^ (2.13) 

with Ce — 1.27 and ~ 0.92 [11]. The damping terms DPi correspond to the 
decohering scatterings, which tend to destroy the coherence of the oscillation. In 
density matrix language such terms appear as relaxation terms suppressing the off- 
diagonals of the density matrix, or equivalently the Pa;,y-components of the Bloch 
vectors. The magnitude of the decohering terms is just half of the corresponding 
scattering rate 

D = lr. (2,14) 

The equation of motion for anti-neutrinos can be found by substituting 

and pi, = to the above equations (the latter condition is true when neutrinos 

are in chemical equilibrium). 

The repopulation term used above is only an approximation for the correct elastic 
collision integral, and unfortunately it breaks the lepton number conservation. This is 
not a serious problem however, and it can be circumvented by introducing an explicit 
lepton number conserving evolution equation for The appropriate equation, 

introduced in [27], is 

OO 

= ^/ dppW^P, - %p,). (2.15) 

0 

Equation (2.15) can be derived from the original equations in the approximation 
where the L-violation due to the the approximate repopulation term is neglected. 
So, when using (2.15) to evolve L, it is guaranteed that the repopulation term does 
not directly affect the lepton asymmetry generation. The accuracy of the method can 
be monitored by comparing the value of derived from (2.15) and from directly 
integrating the neutrino distribution functions. In practice the agreement was always 
to be good. 

Neutrino and antineutrino sectors are extremely strongly coupled through the 
asymmetry in equations (2.2). In order to compute accurately enough 
despite the numerical round-off errors, it is then necessary to avoid computing 
through differences of large quantities, such as naively appear in the equation (2.15). 
To this end we will dehne the “small” and “large” linear combinations of dynamical 
variables in particle and antiparticle sector 

Pt = Pi^ h (2.16) 

For convenience we also separate active and sterile sectors by dehning 

+ Pf = 2ptJk (2.17) 

P± = P„± - P± = 2p±//„. (2,18) 
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In terms of variables (2.17-2.18) the equations take the form 


nq* + r [2/±//„ - F±], ( 2 . 19 ) 

d,Pf = -VP*. (2.20) 

d,Pi = -(H + n)P* - VlPJ - DPt, (2.21) 

d,P±= (V„+V,)Pi + V,P^-\v^(P^-Pf)-DP}, (2.22) 


where we assumed F = F and dehned 


feq = feqiP: h) ± feqiP, “h) (2-23) 

Because 14 = 14, the equation for the asymmetry now reads simply 

OO 

^ I dp pV,p-. (2.24) 

0 

Finally, due to the particular form of the differential operator dt = dt — Hpdp, it 
is possible to reduce our set of partial differential equations to ordinary differential 
equations by changing variables p —> a; = p/T, and t ^ T, and correspondingly 
using the time-temperature relation 

V = -HT. (2.25) 

dt ^ ^ 

The differential operator then becomes simply 

d, ^ (2.26) 

Let us stress that our introducing the small and large linear combinations is an 
absolutely necessary step for obtaining reliable solutions to QKE’s (2.19-2.22,2.24). 
This can best be appreciated from the hgure (3) below, where we show the energy 
spectrum in neutrino and antineutrino sectors for a representative set of oscillation 
parameters; the asymmetry corresponds to the integral over the difference (not even 
visible!) of the neutrino and antineutrino distributions displayed. 


3. Discretization of the momentum variable 

While necessary, introducing parametrization (2.16) is unfortunately not sufficient to 
tackle the momentum dependent QKE’s. Additional difficulties arise due to struc¬ 
tures in the momentum direction (in variable x) and in particular the very narrow 
ones introduced by the resonances. In this paragraph we explain step by step how 
we discretize the momentum variable such that these difficulties can be overcome. 
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First, the momentum variable ranges over the semi-inhnite range x G [0, cx)], 
but in practice we will introduce cut-offs to both ends of spectrum. First, the ultra- 
relativistic neutrino approximation breaks down when x ^ 0, formally appearing as 
a singularity in at a; = 0 in our equations. Second, it is useful to cut out very large x, 
in the exponentially suppressed tail of the distribution. So, we restrict a; to a range 
X G [xminiXmax]) where in practice it suffices to use^ Xmin = 10“^ and Xmax = 100. 

Second, the variable x is not ideal for discretization, because most of the variation 
in the spectrum occurs at small a; < 3. We therefore have mapped the variable x to 
u{x) belonging to range [ 0 , 1 ] by a transformation 


u{x) 


a:(l -f 62 ) - a;extei 

X -f a;ext 


(3.1) 


where ei = ki{l + € 2 ) and 62 = (1 -l- ki)/{k 2 — ki) with ki = a;niin/a^ext and k 2 = 
a^max/a^ext- Apart from a small correction due to the cut-off parameters (3.1) is 
designed so that the extremum point of the thermal momentum distribution Xext 
gets mapped to u ~ 1/2. We will also need the inverse of this function, which is 


x{u) 


Xe^tiu + ei) 
1 -|- €2 — u 


(3.2) 


While clearly a signihcant improvement for binning the momentum variable, (3.1) is 
not clever enough to solve the numerical problems arising from the the pronounced 
structures around the resonant momenta. Fortunately, the resonance positions 
and Xr 2 can be solved analytically from the conditions W = 0 (neutrinos) and 14 = 0 
(antineutrinos): 


+ X — //—resonance 

Xr 2 = + X + //—resonance 


(3.3) 

(3.4) 


where in the case of and Vr neutrinos 


^ = 

X = 




6.60 X 


14C(4)T2 

C(3)M||(lm^| cos 26^0 21 rr 

- — -1.64 X 10® (5m^ Lv 2 cos 26*oT 

14vl2C(4)GFiV.,r3 I i.v 


(3.5) 

(3.6) 


Using (3.1) one easily hnds the corresponding resonant values in the u-variable. 

The resonance widths can be estimated from the distance over which the matter 
mixing angle drops to a some fraction of its maximum value of unity. In this way one 
hnds that the relative resonance width is proportional to the vacuum mixing angle: 

/\ X 

-~ tan 26^0. (3.7) 

_ X 

^Note that even for x = 10“^ still p > 1 throughout our computation, so that the ultrarel- 
ativistic approximation always remains valid. 
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From (3.7) it is clear that for small mixing angles the structures near resonances 
can only be resolved by an extremely hnely spaced momentum grid. For example^ 
with sin^ 26^0 = 10“®, equation (3.7) predicts structures in the scale Ax/x ~ 10“^. 
Assuming then that at least 50 points are needed to describe the resonant region 
accurately, one sees that a linear binning of the x-variable with Xmax = 100 would 
require 5 x 10® grid points! Solving such problem is clearly beyond capabilities of 
any computer we have today. 

Our mapping (3.1) helps the situation by roughly a factor of hundred, but the 
remaining problem would still be way too large. We therefore need to hnd a second 
transformation u u{v) that would redistribute the points such that more grid 
points are clustered around the resonances in the physical momentum variable. There 
are obviously many ways to implement such a transformation. We found it most 
convenient to construct it by using polynomials of the form a + b{v — VrY. Because 
we have two resonance points the mapping function is a little complicated: 

u{v) = av + {a + b(y — Vr^Y) 0{vc — n) + (c + d{v — Vr^)^) 9{y — Uc), (3.8) 

where the two separate mappings around resonant values and ^^2 are joined 
smoothly at the geometric mean of the resonances Vc = \{vr^ -\-Vr 2 )- The linear term 
av was added in order to have a control how steep is the distribution function of the 
grid points. 

The mapping (3.8) has six free parameters a, b, c, d, and Vr 2 , which will be 
determined from boundary conditions at v = 0 and n = 1 and continuity conditions 
at n = Vc- First, the constants a, b, c and d are hxed by the boundary conditions 


m(0) = 0, 

^(^ri) 

u{Vr2) = Ur2-, 

m(1) = 1, 

which results in 


(3.9) 

(3.10) 

(3.11) 

(3.12) 


b = 


d = 


u 


ri 

^3 

ri 


a 

v'i 

ri 


U 


r2 


a 


(1 


-'r2j 


(1 


\2- 


(3.13) 

(3.14) 

(3.15) 

(3.16) 


''T2 > 


Note that are not unknowns, but can be solved analytically from the resonance 
conditions (3.3-3.4) and the inverse mapping (3.2). Finally, the continuity conditions 


lim 

V^Vc + 


U[V] 


= lim u{v) 

V^Vc — 


(3.17) 


®One is driven to such small values of vacuum mixing to find acceptable phenomenology for the 
asymmetry growth [15,32]. 




Figure 1: The mapping u = u{v) is shown with = 4 x 10 T = 11 MeV (solid line) 
and with = 10“^®, T = 15 MeV (dashed line). Mixing parameters are = —0.1 eV^ 
and sin^ 26q = 10“®. 


lim d^uiv) = lim dyu{v). (3.18) 

V^Vc+ V^Vc — 

provide implicit equations for the remaining two unknowns Vn'- 

a + ^b{Vr 2 - Vr^)^ = C+ ^d{Vr^ - (3.19) 

O O 

b = d. (3.20) 

Unfortunately these are sixth order algebraic equations and can only be solved by 
iterative numerical algorithms. 

The novel feature of the above change of variables is that it automatically follows 
the resonances, placing majority of the grid points dynamically in their vicinity. 
As an alternative to our approach, one could think of using library routines with 
automatic grids to take care of this distribution. However, the resonance positions 
evolve as a function of expansion scale of the universe and in particular as a function 
of High point densities are thus needed at different momenta at different times, 
and this movement can be fast. None of the library routines we tried came close to 
being able to cope with these requirements and failed completely to solve the problem. 
With our parametrization on the other hand, the numerical program can be written 
for a hxed grid in the variable v, and the redistribution of points is provided by the 
dynamical equations themselves. 

In Fig. (1) we show the transformations u —> n when jg almost zero (dotted 
line), and when has grown to a snbstantial value (solid line). In the former 
case one cannot distinguish the two resonances by eye, whereas in the second case 
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Figure 2: Density of points (defined as the inverse of the step size) in the physical mo¬ 
mentum variable x for the transforms shown in figure (1). 

the structures have become very distinct. Fig. (2) shows the density of points as a 
function of x (dehned as the inverse of separation of grid points in the x-variable). 


4. Changes in the QKE’s due to the parametrization 


The changes of variables (3.1) and (3.8) affect the partial derivatives appearing in 
Eqns. (2.19-2.22,2.24). These equations have been written assuming the evolution 
of p and L is along constant x contours. However, we have now hxed our grid in 
the variable n, and hence need to find the evolution equations with T and v as the 
independent variables. These are easily found from the old ones by the chain rule 
partial differentiation 


dp{T,x{T,v)) \ ^ /^\ 
dT )v KdTJx 



dx \ 

(—) 



&T)v 

\dxJ T 



dx \ 1 

^du\ 1 


f^\ 


KdxJ 

•.du) T 

Kdv) T 

du\ , 

^dv\ 

f—) 


dTJv^ 

kOuJ T 

\dv ) T 



In the last equality we used the fact that u = u{x) only. The differentials {dTp)x are 
of course just given by the l.h.s. of equations (2.22), and the last coefficient {dvp)T 
is computed numerically as a part of the system of the partial differential equations. 
To be precise, we are using the central derivative formula 


/ dp\ p{v + h) — p{v — h) 

) rp 2h 
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Of the remaining differentials {dvu)T is also easily solved from the parametrization 
Eq. (3.8: 


^ a + 3[b{v — Vri)'^ 0(Vc — v) + d{v — Vr 2 y 0{v — Vc)] 

However, solving the remaining differential {Otu)^ is somewhat more complicated. 
Using equation (3.8) one can write 


(Otu)^ = [dra + dTb{v — — 36(u — d{vc — v) 

+ [OtC + dTd{v — ~ 3(i(u — Ur2)^'9TUr2] — Vc), (4.4) 

where 


dra = drUri — 

drb = ^ (dTUri + (2a- ^^')dTVri 
vf, \ V Vr, j 

OtC = dTUr2 — adTVr2 
drd = —— - — (dTUr2 + ( 2 a — 3-- -^dxVr^ . 

{1-Vr2r\ ^ 1-Vr2'' J 


(4.5) 

(4.6) 

(4.7) 

(4.8) 


Here the only unknowns are the temperature derivatives of the resonances in the 
u-variable Otv^- These can be solved from our original matching conditions (3.19) 
and (3.20), which can formally be written as 


/i(u,,(T),u,,(T),T) =0 (4.9) 

f2{Vr,{T),Vr2{T),T)=0. (4.10) 


Differentiating equations (4.9-4.10) with respect to T results in a linear set of equa¬ 
tions for drVri, which after a little algebra can be brought to the form 

drVr^ = UijdTUrj, (4.11) 

where the matrix Uij depends both on Ur^ and Note that despite the apparent 
complexity of our parametrization, we have been able to reduce everything except 
solving for the resonance parameters to a simple linear algebra. Solving from 
the matching equations at each time step would be very time consuming however. 
This is so in particular because we need to solve them with very high accuracy in 
order for not to lose the continuity of the matching at v = Vc and to avoid inducing 
random numerical errors into the equations. Fortunately, the equation (4.11) itself 
offers the way out of this problem. Namely, we only need to solve from the 
matching conditions (3.19) and (3.20) once, in the beginning of the iteration, as 
their value can later be evolvedhy equations (4.11) as the solution proceeds. In this 
way even solving for the parametrization from v to x becomes part of the dynamical 
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p/T 


Figure 3: Example of neutrino {paa) and antineutrino {paa) spectra as a function oipjT. 
Lines fall on top of each others. Gray line marks the resonance momenta. 

equations. The accumulation of numerical error can be easily traced by solving the 
matching conditions independently for some temperatures. In practice the errors do 
not accumulate at all and it turns out that using the differential equation is by far 
the superior method for solving in comparison to using the matching conditions 
at each temperature step. 

5. Numerical Results 

We have solved the evolution equations numerically for a range of representative 
parameter sets. In hgure (3) we show the neutrino and antineutrino distribution 
functions in a relatively early stage in the evolution for 5rn? = —10 eV^ and sin^ 2^0 = 
IQ-e.i distributions show a very nice thermal shape, and no visible difference 

between particle and antiparticle sectors can be seen. This is so despite the fact that 
at the temperature at which the distributions were plotted, the momenta p ~ 1.8T 
are currently resonant (indicated by the gray line in the plot); the amplitude of 
resonant structure is simply too small to be seen in the scale of the “large variables”. 

The resonance does alter the asymmetry spectrum however. In hgure (4) we show 
the distribution of the difference of neutrino and antineutrino densities p^a — Paa in 
the physical variable x. The prominent and extremely narrow resonance structure 
(one cannot separate the neutrino and antineutrino resonances here) is clearly visible. 
The same spectrum is shown in the integration variable v in hgure (5). Now the 
sharp kink-like resonant structure has broadened into a broad wave, which is easily 
represented by the linearly discretized function in variable v (grid points are shown 
by the open circles). In this example we used just 400 grid points, but essentially 
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Figure 4: Difference of neutrino and antineutrino densities paa — Paa as a function of p/T 
for the same parameters as in figure 3. 
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Figure 5: Same as figure (4), but shown against variable v. 

identical results were found with only 200 points. Moreover, comparing the scales in 
figures (3-4), one sees that the scale of variation of L{x) is ten orders of magnitude 
below unity. So, in order to follow the evolution of accurately when using the 
original variables, one should be able to integrate over distributions with an accuracy 
much better than ten digits! This seems totally impossible, and clearly indicates that 
our use of variables is essential. 

These examples prove the success of the key elements in our approach, namely 
that by separating the small and large variables, and by introducing the novel dynam- 
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ical discretization of the momentum grid, we can model the evolution of the ensemble 
of oscillating neutrinos accurately. Moreover, we can do so using small enough lat¬ 
tices such that the code can be run even in a powerful modern workstation (although 
a single run may then take as long as a couple of days). 

5.1 Regulating the sterile neutrino spectrum 

As one lets the system evolve further the resonance x moves to larger values while the 
integrated decreases. This continues until one reaches roughly the point when 
the main bulk of the active neutrino distribution is getting through the resonance. At 
this point one expects that the total asymmetry begins to oscillate or grow, depend¬ 
ing on the parameters. In our initial runs we did not get that far however, but the 
program crashed due to an accumulation of structures in sterile neutrino spectrum 
with momenta less than XresT. How these structures arise is easy to understand. 
Indeed, the active sector distributions are kept smooth away from the resonances 
by diffusion induced by collisions with the background particles. Sterile states in¬ 
teract only through the effective interaction due to mixing however, and since this 
interaction is neglibly small away from the resonances, the diffusion does not work 
on sterile sector. As a result the moving resonance leaves the sterile neutrino mo¬ 
mentum distribution quite oscillatory for x < Xres- While these developments in the 
sterile sector have little effect on variables in the active sector (because active and 
sterile states are essentially decoupled away from resonances), they still pose a nu¬ 
merical problem: as the resonance moves forward in x, the density of points reduces 
drastically in the region x Xresi until the point that the grid becomes too sparse 
to model the oscillatory structures in the unsmoothed sterile neutrino distribution. 
Eventually the amplitudes of these structures overflow causing numerical instability 
and the solution breaks down. 

The above explanation of the instability also suggests an obvious cure for the 
problem. If one is not interested in studying the sterile neutrino asymmetry spec¬ 
trum, all we need to do is to add a small regulatory interaction term to sterile neutrino 
sector, which will smooth their distribution away from resonances, and yet does not 
alter the evolution of the variables in the active sector. The latter requirement can 
be met by a suitable form of the interaction term and by making it small enough so 
that it does not change the physics in the resonance region. The challenge is to do 
this so that the regulatory term remains efficient enough to numerically stabilize the 
system. To this end we have added the following “sterile repopulation terms” onto 
equation (2.20): 

Rts = rsT {UuJeqix, /i) - pi) (5.1) 

(pii - p- - p-) . (5.2) 

Here p™* is the initial density matrix in the active sector. We used the active sector 
interaction rate T above, and let the parameter control the size of the regulatory 
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Figure 6: The stability of asymmetry evolution against variation of parameter We 
have used values Vg = 0.0001 (solid line), 0.001 (dashed line), 0.01 (dotted line) and 0.1 
(dash-dotted line). Solutions for values r* = 10“^ and r* = 10“^ (long dashed line) are 
already indistinguishable in the figure. Oscillation parameters are = —0.1 eV^ and 
sin^ 200 = 10 “®'^. 

term. The regulators (5.1-5.2) differ from the repopulation terms in the active sector 
in (2.19) in that we have added normalization factors to make sure that the integrals 
over always vanish. (In this way we avoid spurious sterile state equilibration 
even with large control parameter r^.) 

We obviously need to show that there exists a range of values for which the 
results converge and become independent of r*, while the regulator still continues 
to stabilize the solution. That we do reach this goal is evidenced in Fig. (6), where 
we show the evolution of the total asymmetry for = 0 . 1 , 0 . 01 , 0 . 001 , 10 “^ and 
10 “^, for 6nR = — 0.1 eV^ and sin^20o = 10 “®'^. Unsurprisingly, for large there 
is a signihcant effect, but for small the results do converge to the point that no 
visible difference can be seen between the two curves corresponding to two smallest 
r^’s. Let us stress that the mixing parameters were here deliberately chosen such as 
to show a particularly large sensitivity for the regulator. For many other cases that 
we studied it would have been difficult to show any effect at all; for example the 
results displayed in hgure (8) showed no appreciable dependence on the regulator up 
to Vg = 10 . 

We have also studied the stability of our results against changes in various nu¬ 
merical error handling parameters in our code. The stepwise error tolerance for 
example can be relaxed by an order of magnitude from what we actually employed, 
before any visible deviation from the solution can be seen over the entire integration 
range. This is true also for our examples corresponding to the oscillatory solutions 
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in or near the chaotic region. Of course, the numerical accuracy would be eventually 
compromised for parameters for which the asymmetry undergoes very large number 
of oscillations before settling to the growing curve, but as in [17], that does not affect 
our main conclusions. Also cut off parameters x^in and Xmax and the tilt a which 
controls the density of points at the resonance can be varied within reasonable limits 
without any observable effects on the results. We therefore are conhdent that our 
observed pattern of asymmetry oscillations is indeed physical, and not of numerical 
origin. 

5.2 Argument for chaoticity in L-growth 

Having introduced a method that can handle the evolution of momentum dependent 
density matrices and in particular the asymmetry, and having proven its numerical 
stability and accuracy, we now pursue the question of the “chaotic” behaviour of 
the hnal sign of the integrated asymmetry. Note that while the use of word chaos 
is customary in this context, we strictly speaking mean only sensitivity of the hnal 
sign of the asymmetry on initial seed asymmetry and on oscillation parameters®. 

In the context of our earlier work using the momentum averaged equations [17, 
19], we explained what was the physical origin for the sign sensitivity. Quite simply, 
changes in the boundary conditions (initial baryon asymmetry and the oscillation 
parameters) change the direction of the motion of the system, and the topography 
of the attractor solutions in the phase space at the onset of resonance. Roughly 
speaking, in the momentum averaged case one can draw an analogy between the 
system and a somewhat sticky ball initially rolling down a single valley hoor that 
later branches to two, with a low maximum in between, at the resonance point. 
Changing baryon asymmetry would then be analogous to changing the speed and 
direction of the ball at the branching point and changing oscillation parameters to 
that of changing the shape of the valleys themselves. Sign sensitivity arises for those 
initial conditions and valley forms for which the ball makes many oscillations between 
the two degenerate valleys before stickiness (damping) eventually wins and makes it 
to settle into either one of them. 

In the momentum dependent case these structures become smeared out, and to 
draw a mechanical analogy, one should rather think, instead of a ball, of a string of 
varying size beads (largest ones in the middle, corresponding to the peak in the ther¬ 
mal distribution) coupled by a very elastic cord, wiggling down the valleys described 
above. While this is obviously a much more complicated system with qualitatively 
new features, one would still expect to see oscillations in the weighted average posi¬ 
tion of beads at least when the largest ones roll past the branching point (in L when 
the maximum of the distribution crosses the resonance). Moreover, if the average 

®In fact the system does exhibit true chaoticity at a certain level, but the related information 
loss is very small in practice [33]. 
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Figure 7: Neutrino asymmetry oscillations are shown near the boundary of ’’chaotic” 
region using same = —10 eV^ and three slightly different mixing angles sin^ 20o = 
10“®'^ (solid line) sin^ 26 q = (dotted line) and sin^ 26 q = 10“®-®^ (dash-dotted 

line). 


position of the string (L) makes many oscillations over the central hill before stick¬ 
iness (damping) wins, one again should expect that the hnal valley that the string 
settles in (sign of L) strongly depends on initial and boundary conditions. While 
this analog is certainly bears only a crude resemblance to the true system, it should 
help convince the reader of the fact that there is nothing mysterious about the sign 
sensitivity in the L-evolution, but that it is rather something to be expected. 

Obviously, all we need to do to prove the existence of chaoticity, is to hnd exam¬ 
ples of almost identical sets of oscillation parameters for which the hnal asymmetry 
does show signihcantly different oscillation pattern, and leads to a different hnal 
sign. We show such a case in hgure (7). All three curves in the hgure correspond to 
5m? = —10 eV^ but have slightly different vacuum mixing angles. The two curves 
which behave almost alike, have mixing angles sin^ 26q = 10“®'^ (solid curve) and 
sin^26*o = 10“®-®®® (dotted curve) respectively. The dash-dotted curve with a very 
different behaviour corresponds to sin^20o = 10“® ®®. So, if such mixing were ever 
observed and the mixing angle was measured to a precision better than a tenth of 
a per cent, the final sign of the lepton asymmetry and its effect on nucleosynthesis 
could be predicted from our computation. However, if the mixing angle was mea¬ 
sured with less than about one per cent accuracy, we would not be able to say what 
the final asymmetry will be, and therefore what is the precise BBN prediction for 
helium abundance. 

The sign sensitivity does not occur for all oscillation parameters, however. In the 
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Figure 8: Asymmetry evolution in the stable region of oscillation parameters. All curves 
corresponds to mass splitting Sm'^ = —0.1 eV^ and the mixing angles are sin^ 20o = 10“® 
(solid curve), sin^ 26q = 10“®'® (dotted curve) and sin^ 29o = 10“®'^ (dashed curve). 


figure (8) we show results from “stable” region of parameters (same mass difference 
as above, but much smaller mixing angle). In this case, while visible differences in 
the evolution become visible when mixing angle is changed by a factor or three, the 
asymmetry never becomes oscillatory, and the final sign is determined by the initial 
sign of the asymmetry. 

Our first example was chosen from close to the edge of the chaotic region, and just 
as making sin^ 29q smaller made system more stable, increasing it has the opposite 
effect. The number of oscillations and also degree of sensitivity of the final sign to the 
oscillation parameters increases very rapidly with sin^26'o. Nevertheless, the degree 
of sensitivity to parameters is weaker and thus the chaotic region smaller than what 
was found in the momentum averaged treatment [17]. Quantitatively our new results 
agree roughly with the boundaries of the chaotic region inferred in [27], but since the 
numerical work is, even with the improvements we have made, very time consuming, 
we will not attempt a precise mapping of the chaotic region here. 

Let us finally to try to understand a little more deeply what precisely drives 
the oscillations of the total asymmetry in the momentum dependent case. A clue 
can be found from Figure (9) which shows the variable P~ in the region where 
the asymmetry oscillates significantly. (At earlier times, when the total L does not 
oscillate, the structure of the P“-spectrum is similar to that of p~^ shown in Fig. 
(5).) The two new peaks seen outside the resonance grow slowly in time, while their 
sign oscillates rapidly. Because these structures are asymmetric around the resonance 
(in the physical variable x), the integral over P^{x) is nonzero and oscillates. As a 
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Figure 9: Lower panel; Variable P~ as a function of v at temperature T = 27.8 MeV for 
oscillation parameters 5vn? = —lOeV^, sin^ 20o = Upper panel: The coherence 

length (dashed line) and ^m/4 (solid line) for the same parameters. Physical momentum 
at the resonance is Xres = 1-88. 

result also L, which is solved from (2.24), becomes oscillatory. 

What causes these structures, is a delicate interplay between the coherence length 
and the matter oscillation length. The former is simply given by the inverse of the 
damping rate f'coh = 1/D and the latter is 

= (5.3) 

a/ 1 - 2xcos26»o + x® 

where /vac = 4:np/\6m^\ and y = 2p\V\/\Sm‘^\ where \V\ is the magnitude of the 
matter contribution to the W in equation (2.4). Near the resonance the matter 
oscillation length is very long, and in particnlar at high temperatures it exceeds the 
coherence length by a large margin. In such a situation the resonance is strongly 
overdamped, and despite the maximal mixing angle oscillations do not have time to 
develop because the state is continuously projected to the active direction. This in 
part explains why the resonant featnres have so modest amplitudes at early times 
(cf. Fig. (4)). 

In the npper panel of the hgure (9) we have plotted the coherence length and 
£m/4, which corresponds to the distance over which hrst rises to maximnm due 
to oscillation. It is clear that the new structures correspond to the case when the 
matter oscillations first time break the overdamping sitnation. Further away from 
the resonance the damping no more can block the oscillations, bnt since the matter 
mixing angle decreases very fast away from the resonance, the amplitude of the 
oscillation features dies off quickly. As the temperature decreases, the new structnres 
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move away from the resonance, and they also become broader (one can show that 
their separation scales like Sx ~ T~^). This is why their amplitude slowly increases in 
time (with decreasing T) causing larger and larger amplitude oscillations in the total 
asymmetry L. The lesson to bring home from these observations is that the dynamics 
of L-growth and oscillations is very dependent on the detailed hne structure of the 
various components of the density matrix. Such hne details on the hand can only 
be studied in a numerical approach using both separated small and large variables 
(P^) and dynamically adjusting momentum grids. 

6. Conclusions 

Active-sterile mixing, while constrained by the present solar and atmospheric data, 
are nevertheless required if one wishes to incorporate the LSND-neutrino anomaly 
into the neutrino models. Moreover, the active-sterile neutrino mixing would have 
a rich phenomenology in the early universe, where it provides a unique theoretical 
challenge in the form of a system for which the quantum effects play an essential role 
in the kinetic equations. 

Here we have studied the phenomenon of neutrino asymmetry growth in the early 
universe as a result of active sterile neutrino oscillations. Our main new contribution 
is the introduction of a novel method of discretizing the momentum variable such 
that the sharply pronounced structures near the resonances, which are here shown 
to drive the oscillation phenomena, can be treated numerically accurately. The only 
concession to rigor we made along the way to solution of the problem, was adding 
collision terms to regulate the sterile neutrino spectrum away from the resonances. 
However, we demonstrated that our regulators had no effect on the active neutrino 
evolution, and hence for the results presented here. Nevertheless, our treatment 
is obviously not adequate if the precise form of the sterile neutrino spectrum is 
important. This might be the case for the models where sterile neutrinos are invoked 
to provide a non-thermal component for the dark matter [34]. 

We have demonstrated that even tiny changes in the oscillation parameters may 
drastically alter the oscillation pattern, and even change the sign of the hnal asym¬ 
metry. This behaviour does not occur for all oscillation parameters, but instead 
the dependence of the sign of L on oscillation parameters is weak in the “stable” 
region (in general small sin^ 2^0 and large negative 6m^) and strong in the chaotic 
region (in general large sin^ 26^0 and small negative 5w?). These results are in qual¬ 
itative agreement with our earlier hndings, which were based on a momentum av¬ 
eraged treatment [17]. Obviously, if active-sterile mixing were to be observed with 
parameters residing in the chaotic region, it would be, due to inevitable errors in 
the experimentally measured parameters, impossible to reliably determine the sign 
of the neutrino asymmetry created by that mixing in the early universe. Moreover, 
since large neutrino asymmetries (and electron neutrino asymmetry in particular) di- 
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rectly affect, in a si 5 'n(L)-dependent way, the weak interaction rates that determine 
proton-to-neutron ratio in the early universe, this indeterminacy would undermine 
our ability to accurately compute the BBN prediction for light element abundances 
for chaotic oscillation parameters. 

These conclusions hold in the spatially homogeneous calculation. However, the 
sensitivity of L-growth on initial conditions also lends to a speculation that in reality 
the initially inhomogeneous seed asymmetry might give rise to an inhomogeneous 
texture of domains of very large lepton asymmetries with oscillating sign. This 
phenomenon has been shown to occur in a one-dimensional, momentum averaged 
model [35], and could plausibly occur in a realistic three dimensional world. The 
implications of such a scenario for SBBN would obviously be very different, including 
the possibility of very efficient equilibration of sterile neutrinos via MSW-effect within 
the domain boundaries [35,36]. If this was the case, then the asymmetry growth 
mechanism would lead to even stronger bounds on mixing than what is displayed in 
equations (1.2). It is beyond the scope of this work (and the reach of the present 
computers) to study this phenomenon quantitatively in the momentum dependent 
case, however. 

Let us hnally comment on the effects of large homogeneous lepton asymmetries 
for cosmic microwave background radiation (CMBR). The possibility of measuring 
L by the Planck satellite data has been considered for example in [37] and in [38]. 
However, the only effect of L on CMBR comes through the associated fluctuations in 
the energy density, and correspondingly only the total asymmetry has relevance. In 
the oscillation scenarios, such as discussed in this paper, the total asymmetry (here 
the sum of the sterile and active sector asymmetries) is conserved however, and hence 
the oscillation-induced asymmetries, in contrast to the situation with SBBN, would 
have no direct effect on CMBR. 
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